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Abstract. SIBYLL 2.1 is an event generator for hadron interactions at the highest energies. It is commonly 
used to analyze and interpret extensive air shower measurements. In light of the first detection of PeV neutrinos 
by the IceCube collaboration the inclusive fluxes of muons and neutrinos in the atmosphere have become very 
important. Predicting these fluxes requires understanding of the hadronic production of charmed particles since 
these contribute significantly to the fluxes at high energy through their prompt decay. We will present an 
updated version of SIBYLL that has been tuned to describe LHC data and extended to include the production 
of charmed hadrons. 


1 Introduction 

SIBYLL □ is a hadronic interaction model that is widely 
used in air shower simulations. It is available as one of 
the standard hadronic interaction models for high energy 
in the simulation packages AIRES, CORSIKA, CONEX 
and SENECA. SIBYLL is also used for calculating atmo¬ 
spheric lepton fluxes ||2|. 

The cuiTent version of the model is SIBYLL 2.1 El. 
It is designed to allow simulation of hadronic interactions 
in the energy range from ^/s =» lOGeV up to 400 TeV. 
At the time of tuning the parameters of this model, TeVa- 
tron data (^[s ~ 2 TeV) were the highest energy mea¬ 
surements available. In this work we present a new ver¬ 
sion of SIBYLL tuned to accelerator data including those 
from LHC. In addition, this version has been extended to 
include a phenomenological model of the production of 
charmed hardens. 

In the first section we describe the updates of the model 
motivated by LHC data. This includes the refit of the 
cross section parameters, the extension of the fragmen¬ 
tation model to increase baryon pair production, and the 
update of the parton distribution functions. In the second 
section we describe the model of charm production, how 
the parameters are adjusted to describe data and what con¬ 
clusions can be drawn from applying the model in calcu¬ 
lations of the inclusive flux of atmospheric leptons. 

2 LHC updates 

Before discussing the update of the model it is worthwhile 
to mention that SIBYLL 2.1 already describes the general 
characteristics of hadronic interactions at 7 TeV remark¬ 
ably well (see dashed blue histogram in Eig. |^or the re¬ 
view by d’Enterria et al. PI). 



Figure 1. Inelastic p-p cross section in SIBYLL. The updated 
cross section is shown in blue, the old version is in black. The 
red squares are the measurements by TOTEM Q. The black 
diamond at the highest energy is the estimate from the Auger 
Observatory HIH. The second energy axis shows the equivalent 
laboratory energy for p-p interactions as applicable to air shower 
detectors (one proton at rest). The measurement ranges of the 
IceTop air shower array (S) and the Pierre Auger Observatory Qj 
are indicated by black lines. 


2.1 Proton proton cross section 

The hadron-proton cross section in SIBYLL follows from 
unitarizing hard and soft cross section contributions, sep¬ 
arated by an energy dependent cutoff in and terms 
due to diffraction dissociation. More details on the struc¬ 
ture of the model can be found in the publication for ver¬ 
sion 2.1 lO. 
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Measurements at LHC suggest (see Fig. that 
SIBYLL 2.1 overestimates the cross section at high ener¬ 
gies. The inelastic cross section measured in the TOTEM 
experiment, which has the highest precision for a measure¬ 
ment of the total cross section at LHC, is 73.5;^[ 4 mb fSl 
whereas SIBYLL predicts 80 mb. The rise of the pp cross 
section beyond 1 TeV is mainly driven by hard parton scat¬ 
tering (hard minijets). 

In SIBYLL an eikonal approximation is used to com¬ 
bine the parametrization of soft scatterings with the pertur¬ 
bative calculation of the minijets into an unitary amplitude, 
which then defines the total and elastic cross sections. The 
size of the soft and hard contributions in this formalism 
depends on the size of the particular cross section and the 
prohle function. 

In order to make the inelastic cross section compatible 
with the TOTEM result without changing the hard cross 
section (calculated within QCD), the profile function of 
the distribution of the hard partons in transverse (impact 
parameter) space has been made more narrow so that pe¬ 
ripheral collisions are less likely to produce minijets. 

The downside of this approach is that central colli¬ 
sions now exhibit very high densities of interacting par- 
tons (profile functions are normalized), which means that 
some events will have a large number of minijets and 
consequently a large number of final state particles pro¬ 
duced. This effect will produce a tail in the multiplicity 
distribution that is not observed in data. However cen¬ 
tral collisions are rare so the average multiplicity and most 
other observables are still compatible with the measure¬ 
ments n. 

Since our goal is a model capable of describing inter¬ 
actions a decade and more higher in center-of-mass en¬ 
ergy, the effects of high parton densities have to be con¬ 
sidered, even if the mean multiplicity still agrees with cur¬ 
rent experiments. A microscopic model of parton density 
saturation could limit the number of scatterings in central 
collisions and thereby repair the multiplicity. In the cur¬ 
rent model, saturation is implemented only in an impact 
parameter independent way as an energy-dependent lower 
Px-cutoff for the minijets. 

In addition to changing the hard profile function we 
adjust the parameters of the soft cross section parametriza¬ 
tion to ht the p-p and p-p cross section data. 

Since the proton prohle also enters the meson-nucleon 
cross sections, we reht the parametrization of the soft con¬ 
tribution there as well. 

The resulting cross section is shown in Eig.[2as a blue 
line. The old cross section for comparison is shown as a 
black solid line. The data point of the highest energy is the 
estimation of the p-p cross section from air shower mea¬ 
surements at the Auger Observatory at energies of about 
■\fs - 57 TeV Elll. This value has not been used to ht 
the cross section in SIBYLL and therefore can be seen as 
an indication that the extrapolation by the model is reason¬ 
able. 


laboratory energy (eV) 

10“ 10“ io“ lo” lo” 10“ io“ 



Figure 2. Average antiproton multiplicity as a function of center- 
of-mass energy. The low energy data are a compilation of hxed 
target and ISR experiments that cover the full phase space or 
were extrapolated to full phase space |9(. The CMS data HOI are 
taken in a phase space region with \y\ < 1.0. PHENIX im data 
are taken in the range I 77 I < 0.35. The prediction by the models 
are shown for the full and CMS phase spaces only. SIBYLL 2.1 
is shown as dashed line, the updated version as solid line. 


2.2 Baryon production 

Particle production in interaction models primarily de¬ 
pends on the implementation of the fragmentation process. 
Eragmentation is a non-perturbative process so the rates of 
particle production cannot be calculated from hrst princi¬ 
ples, which means the parameters in the model have to be 
set by comparison with experiment. 

In SIBYLL string fragmentation m is used as 
the fragmentation model. The string model simplifies 
hadronization by assuming a uniform energy density in 
the color held stretched between two partons which even¬ 
tually is split in two by quark-antiquark pair production. 
The splitting is continued until the remaining energy is just 
enough to form two hadrons. Baryons are produced by 
introducing diquark-antidiquark pairs instead of qq pairs. 
The probability of producing a diquark pair rather than a 
quark pair (Fq/qq) in a string breakup is the parameter that 
controls baryon pair production. In version 2.1 it is set to 
0.04. 

Eor simplicity, only two string classes are distin¬ 
guished in SIBYLL: the 2 string conhguration for the 
2 —> 2 sea parton scattering and two single strings con¬ 
necting valence quarks/diquarks. The essential difference 
between the two is that the latter conhguration has valence 
Havor attached to the string ends, where as the former is 
in total havor neutral. This distinction is necessary to de¬ 
scribe the differences between the forward/backward re¬ 
gions and the central region of phase space. 

The result of this treatment of baryon production in 
SIBYLL 2.1 for the antiproton multiplicity is shown in 
Fig. 1^ as dashed black lines together with a compilation 
of data. The multiplicity for full phase space, typically 
measured in hxed target experiments at low energies, is 
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shown in the upper set of lines whereas the multiplicity 
in the central region (It/I < 2), the region typically acces¬ 
sible in collider experiments, e.g. CMS cni, is shown in 
the lower set. The current model describes the threshold 
at low energies well but is not capable of describing the 
central, high energy data at the same time. 

In order to allow for a meaningful extrapolation to 
high energies, instead of introducing an arbitrary energy 
dependent parametrization for f^q/qq, one can relate the 
baryon production frequency to soft and semihard inter¬ 
actions, whose energy dependence is different, to increase 
the baryon production mainly at high energy. With the 
minijet cross section being derived from perturbative QCD 
the extrapolation to higher energy is then given by the 
model and, at low energy, threshold effects due to the large 
mass of the baryon pairs are important. 

Furthermore minijets mostly produce particles in the 
central region which is exactly where the high energy data 
by CMS reveal a deficit for SIBYLL 2.1. This assumption 
is supported by the observation of the ratio of antiprotons 
to charged pions compared to the central charged multi¬ 
plicity (see e.g. Fig. 15 in Ref IfTOl I. 

The simplest possible coupling of the diquark produc¬ 
tion parameter to minijets is to choose a different but fixed 
value of Fq/diq in the fragmentation of minijets in compar¬ 
ison to all other fragmentation processes. The resulting 
model describes the data much better (solid blue line in 
Fig .|^, especially in the central region. 

Measurements of baryon production at LHC energies 
that cover the forward phase space could test the assump¬ 
tions made in this model. 


2.3 Transverse momentum of minijets 

In SIBYLL 2.1 the momentum fractions that determine the 
kinematics of the minijets are taken from an effective par- 
ton density function ifTSll 

fix) = gix) + ^ [qix) -I- qix)] , (1) 

where gix) and ^(x) are parametrized according to Eichten 
et al. (EHLQ) lfT4l and the quark distribution function in¬ 
cludes contributions from three light flavors (m, d, and s) 
and the valence quarks. 

In the updated version the same effective parton distri¬ 
bution function (PDE) is used but the quark and gluon con¬ 
tributions are sampled from the same PDE parametriza- 
tions (GRV ifTTlfTSll I that are used in the calculation of the 
hard minijet cross section. 

The main difference between these parametrizations 
is the behavior at low x which, in the case of the GRV 
parametrization, is much steeper. 

In combination with the correction of a mistake in the 
definition of the /t™" the steeper PDEs give a better de¬ 
scription of the spectra in the range of intermediate trans¬ 
verse momenta (2 - 5GeV/c) than in SIBYLL 2.1, see 
Fig.[3l 


- Sibyll 2.3rcl - Sibyll 2.1 



Figure 3. Inclusive cross section for charged particles as function 
of the transverse momentum. The results obtained with the old 
and new versions of SIBYLL are compared with CMS data at 
different c.m. energies QSIIIII. 


2.4 Other updates 

Other general and more technical aspects of the model 
that have been updated but are not discussed here are: the 
transverse momentum acquired in the soft scattering of va¬ 
lence quarks as well as in the string fragmentation is now 
sampled from an exponential transverse mass distribution 
rather than a Gaussian as in the previous version. 

Another aspect of direct importance to air shower pre¬ 
dictions is the enhanced forward production of vector 
mesons with respect to pseudoscalar mesons (pions) in 
meson nucleon interactions HD. Since this mechanism 
has a large influence on muon production in air showers it 
has been implemented in the new version of SIBYLL. 

Eurthermore the implemented Glauber model for the 
calculation of the different cross sections (total, elastic, 
diffractive, and quasi-elastic) in hadron-nucleus interac¬ 
tions has been extended to include a consistent treatment 
intermediate low-mass excitations, leading to enhanced 
screening effects ll20l . 

2.5 Comparison to data 

To show the compatibility of the updated model with ex¬ 
perimental data we look at the charged particle pseudora¬ 
pidity distribution. The advantage of this observable is that 
it is very sensitive to the details of the parton level inter¬ 
action structure and kinematics (ujets, JC;) as well as to the 
subsequent fragmentation process (dV^^j^^/dr/). 

The changes introduced in the cross section are ex¬ 
pected to increase the central multiplicity at energies be¬ 
yond 1 TeV whereas the increased baryon production, due 
to the higher mass of baryonic particles, can be expected 
to lead to an overall decrease in the multiplicity. Eig. 
shows that both effects approximately cancel one another 
at LHC energy and that also the updated model (solid blue 
histogram) describes the CMS data well. 
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Figure 4. Pseudorapidity distribution of charged particles. The 
data are from NA22 |2D,UA5 G2l,CDF (23] and CMS QS). The 
prediction by SIBYLL 2.1 is shown by the dashed line, the one 
for the updated model by the solid line. 


3 Charm quark extension 

SIBYLL 2.1 is limited to the production of particles con¬ 
taining u, d, and s quarks. In version 2.2c of SIBYLL 
it was shown that a simple phenomenological extension 
of the fragmentation model, based on the family con¬ 
nection between strange and charmed hadrons, can ac¬ 
count for the production of charmed particles at low en¬ 
ergy (241. In this approach the normalization is set by 
the rate at which charm quarks appear relative to strange 
quarks Pc - 0.004. 

3.1 New charm model 

Due to the high mass of the charm quark the production 
of charmed hadrons in the fragmentation process is sup¬ 
pressed by a large factor. Instead the dominant channel 
is the direct production of charm quarks in parton-parton 
scattering. In the context of QCD the leading contribu¬ 
tion gg cc is often referred to as QCD gluon-gluon 
fusion Il25l . The momentum transfer of the reaction due to 
the charm quark mass > Q^i„ ~ Inic means that the 
process can be expected to be calculable within perturba¬ 
tion theory. 

The SIBYLL event generator includes only the dom¬ 
inant terms of hard parton-parton scattering at high en¬ 
ergy and does not distinguish between the hadronization 
of the different parton configurations. All parton-parton 
scattering processes are fragmented into hadrons through 
an unflavored two string configuration, similar to 2 scat¬ 
tered gluons (usually referred to as hard minijets). 

To account for the dominating hard scattering contri¬ 
bution the charm quark fraction is increased in the frag¬ 
mentation of the hard minijets. To keep the threshold 
behavior at low energy the charm quark fraction is sup¬ 
pressed exponentially in low mass strings. Specifically 



where s is the invariant mass of the scattering partons and 
meff = 20GeV^ is the effective mass scale. To account 
for string configurations of higher order charm production 
is not limited to the end of the strings but extends over the 
whole string. This part of the phenomenological model for 
charm production is referred to as perturbative component. 

Next to the dominant contribution from hard scatter¬ 
ing, experiments have shown that there is an asymmetry in 
charm production in the fragmentation region (i.e. at large 
2 Cf) ESI E3, which suggests a contribution from charm 
production in soft interactions. Two models, which can 
be used to explain this forward production of heavy flavor, 
are the intrinsic charm model 1281 and the^avor excitation 
model 12^ . 

In SIBYLL we chose a model which could represent 
either mechanism by adding charm quarks to any string 
attached to soft scattered partons as well (non-perturbative 
component). These will include valence strings which, 
due to the large momentum fraction and the attached fla¬ 
vor of the valence quarks, are able to produce the observed 
asymmetry at large xp. 

3.2 Tuning the charm parameters 

The values of the parameters in Eq. are adjusted sepa¬ 
rately for the perturbative and non-perturbative contribu¬ 
tion. The perturbative part is tuned to describe the p_p- 
spectra of D mesons measured by the ALICE ll^ and 
LHCb ifSTl experiments in central phase space, since this is 
where its contribution is expected to be dominant (Eig.|^. 
The parameters for the soft contribution are set to account 
for the missing production at low energies (Eig.[7|l. 

In Eig.j^the cross section for inclusive charm produc¬ 
tion is shown as a function of the center-of-mass energy. 
The ALICE data include an extrapolation from central to 
total phase space. The cross section for D meson produc¬ 
tion that is measured directly by ALICE is shown by the 
lower blue points and lines. The dotted line represents the 
inclusive D meson production cross section without sub¬ 
tracting the decays of resonances of higher mass, e.g. D*. 
It is shown here because the low energy measurements are 
not corrected for this either. 

The resulting model correctly describes the rise of the 
inclusive charm cross section with energy and reproduces 
the spectra at different energies. 

3.3 Discussion 

Charmed particles were introduced to the model because 
they are expected to contribute significantly to the inclu¬ 
sive flux of atmospheric muons and neutrinos at high en¬ 
ergy The inclusive fluxes can be obtained by solving 
the corresponding cascade equation ll40ll . The results de¬ 
pend on the spectrum weighted moments 



where xp is the energy fraction of the considered final 
state particle in the laboratory frame, ^ is the hadronic 
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Figure 5. Inclusive charm and £)-meson cross sections as a function of c.m. energy. The data at low energy are D-meson cross sections 
in fixed target experiments (261132}®). The measurements at the highest energies are cc from ALICE (30l|35l. Here data are shown 
extrapolated to full phase space (red circles) and visible only (blue empty squares). At intermediate energies the data taken at RHIC by 
the STAR |3^ and PHENIX 1371 experiments are shown (also extrapolated). The model prediction for the inclusive cc cross section is 
shown by the solid line, the prediction for the production of D-mesons is shown by the dotted line. 


- Sibyll 2.3rcl 



Figure 6. Transverse momentum spectrum of different types of 
£)-mesons in the rapidity interval 4.0 < y < 4.5. Data were taken 
at - 7 TeV with the LHCb detector 1^ . 


production spectrum, and y is the power law index of 
the integral all-nucleon spectrum of cosmic rays. With 
1.7 < y < 2.3 it is evident from Eq. |^that the contri¬ 
bution to the lepton flux is largest for charm production at 
large xl- 



Figure 7. Feynman-jc spectra of charged D mesons in p-p fixed 
target interactions with PLab = 400 GeV 1321 and 800 GeV 1331 . 


Unfortunately, particle production at large xl is diffi¬ 
cult to study at high energy. So far there are no experi¬ 
ments that cover this part of phase space and are capable 
of particle identification (PID) at the same time. The most 
forward detector at the LHC with PID capabilities is LHCb 
(2.5 < y < 4.5). For SIBYLL the contribution from this 
phase space to Zpo is only about 10 % at -\As = 7 TeV (in- 
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Figure 8. Weighted spectrum for O-mesons in SIBYLL at 
^/s = 7 TeV. The contributions from the perturbative and non- 
perturbative model components are shown by the blue and red 
lines, respectively. Note the negligible contribution to the energy 
spectrum from the phase space covered by the LHCb experiment 
(2.5 < y < 4.5, green line). 



Figure 9. Comparison of the weighted spectrum of £)-mesons in 
SIBYLL and the MRS model (38]. The MRS spectra are shown 
for different fragmentation assumptions. 


tegrating the green line in Fig. [^. The prediction of the 
contribution of charmed particles to the inclusive neutrino 
and lepton fluxes therefore are not well constrained by the 
measurements at the LHC. In addition, large-XL produc¬ 
tion is dominated by soft, non-perturbative processes, so 
the prediction can not be well constrained by theoretical 
arguments either. 

In Fig. 1^ a comparison of the weighted energy spec¬ 
trum, i.e. the integrand in Eq.[^ for D mesons in p-p inter¬ 
actions between the model by Martin, Ryskin and Stasto 
(MRS) El] and SIBYLL is shown. The MRS model is 
a perturbative calculation of the charm production that is 
extended to low momentum fractions using additional as¬ 
sumptions and accounting for saturation. The energy of 
the comparison is 7 TeV and the index, with which the 
spectrum is weighted, is y = 2. To compare the shapes 


of the distributions the models are scaled such that MRS 
and the perturbative component in SIBYLL are equal at 
xf = 0.19. One can see that the MRS model and the 
perturbative component in SIBYLL show a similar be¬ 
havior. The main difference between the models is that 
SIBYLL predicts additional charm production from the 
non-perturbative component that is dominating in the for¬ 
ward direction. 

A detailed calculation of the atmospheric lepton fluxes 
using the model discussed here can be found in our sec¬ 
ond contribution ATI to this conference. In that paper, the 
role of the all-nucleon spectrum and the atmosphere are 
discussed as well. 

It should be mentioned that the entire discussion here 
was focused on proton-proton interactions. What mat¬ 
ters for the atmospheric fluxes is the charm production 
in nucleon-nucleus interactions. In principle, the model 
is implemented such that central (perturbative) charm pro¬ 
duction should scale approximately with the number of bi¬ 
nary interactions, while forward charm production scales 
with the number of projectile participants. In practice 
these scaling expectations are not really satisfied because 
of additional energy-momentum constraints. Given the 
strong dependence of the atmospheric fluxes on the for¬ 
ward production nuclear screening effects in this region 
could have a large effect. For central production, the mea¬ 
surements confirm the binary scaling 1421 . 

4 Summary and Outlook 

An improved version (2.3rcl) of the hadronic interaction 
model SIBYLL has been presented. The current status of 
the update of the p-p cross section, the extension of the 
fragmentation model to describe increased baryon produc¬ 
tion, and the new charm production model have been de¬ 
scribed in more detail. The perturbative component of the 
charm model was found to be compatible with the analytic 
MRS calculation. It was also shown that the experimental 
data currently available on charm production do not di¬ 
rectly restrict the predictions for the inclusive muon and 
neutrino fluxes. Only indirectly, by comparing model pre¬ 
dictions with charm measurements in phase space regions 
covered also at colliders, constraints on atmospheric lep¬ 
ton fluxes can be derived. 

In the future we plan to estimate how large the uncer¬ 
tainty in the atmospheric fluxes due to the limited phase 
space coverage of the measurements is. This can be 
achieved by looking for a set of parameters in the charm 
model that either minimizes or maximizes the forward 
charm yield while still being compatible with experimen¬ 
tal data. 
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